biomod2建模

library(biomod2);
library(raster);
library(RColorBrewer);
library(dismo);

setwd("~/AnalysisFolder/Localities/")
#Get data points
points <- read.csv(file = "GadMa.csv", header = T);
points <- cbind(points, rep.int(1, length(nrow(points)))); #Adds another column indicating these are presence points
colnames(points) <- c("Species", "X", "Y", "Response");

#Get environmental variables
setwd("~/AnalysisFolder/Ms/GadMa/");
envtList <- list.files(pattern = ".asc");
envt.st <- stack(envtList);

#Get projection variables
setwd("~/AnalysisFolder/EnvironmentalData/Future/rcp8_5/2100/")
projectionList <- list.files(pattern = ".asc");
proj.st <- stack(projectionList);

#Setting up data file for Biomod2
bmData <- BIOMOD_FormatingData(resp.var = points[,4],
                               resp.xy = points[,2:3], 
                               resp.name = as.character(points[1,1]),
                               expl.var = envt.st,
                               PA.nb.rep=1
);

#Setting up Maxent run
myBiomodOption <- Print_Default_ModelingOptions();
myBiomodOption@MAXENT.Phillips$path_to_maxent.jar = paste(system.file(package="dismo"), "/java", sep='');
myBiomodOption@MAXENT.Phillips$memory_allocated = 2048; #Allocates 2048 MB/2 GB of memory to modeling
myBiomodOption@MAXENT.Phillips$maximumiterations = 10000;
myBiomodOption@MAXENT.Phillips$threshold = F;
myBiomodOption@MAXENT.Phillips$hinge = F;
myBiomodOption@MAXENT.Phillips$visible = F;
myBiomodOption@MAXENT.Phillips$beta_lqp = .95;

#Running Maxent
setwd("~/AnalysisFolder/TestModelRun/R/")
myMaxentModel <- BIOMOD_Modeling(data=bmData,
                                    models=c('MAXENT.Phillips'),
                                    models.options=myBiomodOption,
                                    NbRunEval=10,
                                    do.full.models = F,
                                    DataSplit=50,
                                    models.eval.meth = c('KAPPA','TSS','ROC'),
                                    SaveObj = T
);

#Ensemble of all models--combines model runs using a user-selected evaluation metric
myMaxentEnsemble <- BIOMOD_EnsembleModeling( modeling.output = myMaxentModel,
                                   chosen.models = 'all',
                                   em.by = 'all',
                                   eval.metric = c('TSS'),
                                   eval.metric.quality.threshold = NULL,
                                   models.eval.meth = c('TSS','ROC','KAPPA'),
                                   prob.median = TRUE )

#Projecting your model to the present
myBiomodProjPres <- BIOMOD_Projection(modeling.output = myMaxentModel,
                                    new.env = envt.st,
                                    proj.name = 'Present',
                                    selected.models = 'all',
                                    compress = 'gzip',
                                    clamping.mask = T,
                                    output.format = '.grd',=
                                    do.stack=T
);

mod_projPres <- get_predictions(myBiomodProjPres);
presentResult <- calc(mod_projPres,fun = median); #Choose whatever descriptive statistic you'd like
plot(presentResult);
writeRaster(presentResult, filename = "GadusMacrocephalusPresent", format = "ascii", overwrite = T);

#Projecting the ensemble model in the present
myBiomodProjPresEnsemble <- BIOMOD_EnsembleForecasting(myMaxentEnsemble,
                            projection.output = myBiomodProjPres,
                            selected.models = 'all',
                            compress = 'gzip'
);=
mod_projPresEnsemble <- get_predictions(myBiomodProjPresEnsemble);
presentEnsembleResult <- mod_projPresEnsemble[[2]] #This is the median model ensemble
plot(presentEnsembleResult);
writeRaster(presentEnsembleResult, filename = "GadusMacrocephalusPresentEnsemble", format = "ascii", overwrite = T);

#Projecting your model to the future
myBiomodProj2100 <- BIOMOD_Projection(modeling.output = myMaxentModel,
                                    new.env = proj.st,
                                    proj.name = 'In2100',
                                    selected.models = 'all',
                                    compress = 'gzip',
                                    clamping.mask = T,
                                    output.format = '.grd',
                                    do.stack=T
);

mod_proj2100 <- get_predictions(myBiomodProj2100);
result2100 <- calc(mod_proj2100,fun = median); #Choose whatever descriptive statistic you'd like
plot(result2100);
writeRaster(result2100, filename = "GadusMacrocephalus2100", format = "ascii", overwrite = T);

#Projecting the ensemble model in 2100
myBiomodProj2100Ensemble <- BIOMOD_EnsembleForecasting(myMaxentEnsemble,
                                                       projection.output = myBiomodProj2100,
                                                       selected.models = 'all',
                                                       compress = 'gzip'
);
mod_proj2100Ensemble <- get_predictions(myBiomodProj2100Ensemble);
ensembleResult2100 <- mod_proj2100Ensemble[[2]] #This is the median model ensemble
plot(ensembleResult2100);
writeRaster(ensembleResult2100, filename = "GadusMacrocephalus2100Ensemble", format = "ascii", overwrite = T);


#Evaluating models
## Variable response curves
response.plot2(models = BIOMOD_LoadModels(myMaxentModel, models='MAXENT.Phillips'),
               Data = get_formal_data(myMaxentModel,'expl.var'),
               show.variables= get_formal_data(myMaxentModel,'expl.var.names'),
               do.bivariate = FALSE,
               fixed.var.metric = 'median',
               col = brewer.pal(10, "Spectral"),
               legend = TRUE,
               data_species = get_formal_data(myMaxentModel,'resp.var')
);

##Varible contributions; only for Maxent, not possible for other models
forSetup <- read.csv(file = paste("~/AnalysisFolder/TestModelRun/R/GadusMacrocephalus/models/1457123817/GadusMacrocephalus_PA1_RUN1_MAXENT.Phillips_outputs/maxentResults.csv", sep = ""), header = T)#Choose the appropriate model folder with the seed of the analysis you want
variableContributions <- matrix(data = NA, nrow = length(forSetup[, grep('.contribution', names(forSetup))]), ncol = 10);
rownames(variableContributions) <- names(forSetup[, grep('.contribution', names(forSetup))])
colnames(variableContributions) <- c("Run1", "Run2", "Run3", "Run4", "Run5", "Run6", "Run7", "Run8", "Run9", "Run10")
variablePermutationImportance <- matrix(data = NA, nrow = length(forSetup[, grep('.permutation.importance', names(forSetup))]), ncol = 10);
colnames(variablePermutationImportance) <- c("Run1", "Run2", "Run3", "Run4", "Run5", "Run6", "Run7", "Run8", "Run9", "Run10")
count <- 1;
while (count <= 10){
  temporary <- read.csv(file = paste("~/AnalysisFolder/TestModelRun/R/GadusMacrocephalus/models/1457123817/GadusMacrocephalus_PA1_RUN", count, "_MAXENT.Phillips_outputs/maxentResults.csv", sep = ""), header = T);
  variableContributions[,count] <- unlist(unname(temporary[, grep('.contribution', names(temporary))]))
  variablePermutationImportance[,count] <- unlist(unname(temporary[, grep('.permutation.importance', names(temporary))]))
  count <- count + 1;
}
write.csv(variableContributions, "VariableContributions.csv", quote = F);
write.csv(variablePermutationImportance, "VariablePermutationImportance.csv", quote = F);

##Calculate MESS for 2100
mess2100 <- mess(proj.st, rasterToPoints(envt.st)[,-1:-2],-9999);
writeRaster(mess2100, filename = "GadusMacrocephalus2100MESS", format = "ascii", overwrite = T);

##Create dataset for evaluation
###"Cutoff" gives threshold to optimize evaluation metric, "Sensitivity" and "Specificity" are based on this threshold
myMaxentModelEval <- get_evaluations(myMaxentModel, as.data.frame = F); 
write.csv(myMaxentModelEval["TSS",,,,],"TSSEvaluation.csv", quote = F);#Results for True Skill Score
write.csv(myMaxentModelEval["ROC",,,,],"TSSEvaluation.csv", quote = F);#Results for Area Under the Curve
write.csv(myMaxentModelEval["Kappa",,,,],"TSSEvaluation.csv", quote = F);#Results for Cohen's Kappa

biomod2成熟模型

## 成熟的做法：
## 评估物种建模过程中建模的拟合优度：
library(pacman)
p_load(sp,raster,rgeos,rgdal,rasterVis,raster,fs,sf,tidyverse,
       spatstat,maps,maptools,SDMtune,blockCV,sdm,biomod2)

setwd("./xh/")
sa <- read.csv("./four_4kmpoints_xh/xh_sa.csv")[,2:3]
oc <- read.csv("./four_4kmpoints_xh/xh_oc.csv")[,3:4]
na <- read.csv("./four_4kmpoints_xh/xh_na.csv")[,3:4]
as <- read.csv("./four_4kmpoints_xh/xh_as.csv")[,3:4]
### 使用filter-all建模####
loss <- paste0("C:/Users/Administrator/Desktop/xh/f2_envs/",c("BIO4","BIO7","BIO12","BIO11","BIO15","BIO17","GHI"),".tif")
env <- stack(loss)

occall <- rbind(sa,oc,na,as) %>% data.frame()
occall1 <- raster::extract(env,occall) %>% cbind(occall) %>% na.omit()
occ <- occall1[,8:9]
occ

occ_pr <- SpatialPointsDataFrame(coords=occ, data=occ, proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
occ_data <- BIOMOD_FormatingData(
  resp.var = rep(1, nrow(occ_pr)),
  resp.xy = occ,
  expl.var = env,
  resp.name = "occall",
  PA.nb.rep = 1,
  PA.nb.absences = 10000,
  PA.strategy = 'random',
  na.rm = TRUE
)

xh_opt <- 
  BIOMOD_ModelingOptions(
    GLM = list(type = 'quadratic', interaction.level = 1),
    GBM = list(n.trees = 1000),
    ANN = list(),
    RF = list(ntree =1000),
    MAXENT.Phillips = list(hinge=FALSE,threshold=FALSE,
                           memory_allocated = 4096,
                           betamultiplier = 4)
  )

xh_models_global_filter <- 
  BIOMOD_Modeling(
    data = occ_data,
    models = c("GLM", "GBM", "RF","ANN","MAXENT.Phillips"),
    models.options = xh_opt ,
    NbRunEval =5, #NbRunEval    integer, number of Evaluation run.
    DataSplit = 70, ##数据分割，80%
    VarImport = 0,
    models.eval.meth = c('KAPPA', 'TSS', 'ROC'),
    modeling.id = "FILTER_GLOBE_global",
    do.full.models = FALSE
  )

sink("./SDM/FILTER_GLOBE/single_run_global_filter.txt")
get_evaluations(xh_models_global_filter)
sink()

xh_global_filter_ensemble_models <- 
  BIOMOD_EnsembleModeling(
    modeling.output = xh_models_global_filter,
    em.by = 'all',
    eval.metric = 'TSS',
    eval.metric.quality.threshold = 0.6,
    models.eval.meth = c('TSS','ROC','KAPPA'),
    prob.mean = TRUE,
    prob.cv = TRUE, 
    committee.averaging = FALSE,
    prob.mean.weight = TRUE,
    VarImport = 0
  )



##查看建模的评估结果：
## 单物种 建模结果：
pdf("./SDM/FILTER_GLOBE/global_filter_eva.pdf")
FILTER1 <- models_scores_graph(  ##输出多模型比对结果可视化：
  xh_models_global_filter, 
  by = "models", 
  metrics = c("ROC","TSS"), 
  xlim = c(0.5,1), 
  ylim = c(0.5,1)
)

FILTER2 <- models_scores_graph( ##输出run1和run2以及mean的结果；
  xh_models_global_filter, 
  by = "cv_run" , 
  metrics = c("ROC","TSS"), 
  xlim = c(0.5,1), 
  ylim = c(0.5,1)
)

## 查看整合建模的结果：
FILTER3 <- models_scores_graph( ##输出run1和run2以及mean的结果；
  xh_global_filter_ensemble_models, 
  by = "models" , 
  metrics = c("ROC","TSS"), 
  xlim = c(0.5,1), 
  ylim = c(0.5,1)
)
dev.off()
sink("./SDM/FILTER_GLOBE/xh_global_filter_ensemble_models.txt")
get_evaluations(xh_global_filter_ensemble_models)
sink()

# FILTER_GLOBE_models_proj_currente <- BIOMOD_Projection( modeling.output = xh_models_global_filter,
#                                                      new.env = env,
#                                                      proj.name = "current", #projection to analogous climates in occupied range
#                                                      binary.meth = "TSS",
#                                                      output.format = ".img",
#                                                      selected.models='all',
#                                                      build.clamping.mask=FALSE ,
#                                                      do.stack = TRUE )
FILTER_GLOBE_ensemble_models_proj_currentw <- 
  BIOMOD_EnsembleForecasting( EM.output = xh_global_filter_ensemble_models,
                              new.env = env,
                              proj.name= "global_filter_all",
                              do.stack =TRUE )
pdf("./SDM/FILTER_GLOBE/global_filter_ense.pdf")
plot(FILTER_GLOBE_ensemble_models_proj_currentw)
dev.off()
writeRaster(FILTER_GLOBE_ensemble_models_proj_currentw@proj@val[[1]],"./SDM/FILTER_GLOBE/global_ense_mean.tif")
writeRaster(FILTER_GLOBE_ensemble_models_proj_currentw@proj@val[[2]],"./SDM/FILTER_GLOBE/global_ense_cv.tif")
writeRaster(FILTER_GLOBE_ensemble_models_proj_currentw@proj@val[[3]],"./SDM/FILTER_GLOBE/global_ense_wmean.tif")